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ABSTRACT 



Approximation of the differential equation of the single 



degree of freedom, damped oscillator by the first order central finite 
difference equations and two Hermitian finite difference equations is 
investigated. An error analysis is made between the solution to the 
differential equation and the solutions to the finite difference 
approximations for various values of damping. The results of the error 
analysis indicate that the Hermitian approximations have less error 
than the first order difference approximation for the same size of 
increment. Furthermore, the employment of the Hermitian method does 
not materially increase the execution time on the digital computer over 
that required by the first order difference equations. Thus the 
Hermitian equations are superior to the first order finite difference 
equations for the damped oscillator problem. 
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TABLE OF SYMBOLS 



arbitrary constant (in) 

viscous damping coefficient (1^^ sec/in) 

base of the natural system of logarithms (2.71828) 






finite difference station number 

spring constant (Ib^/in) 

2 

mass (lb- sec ) 

^in 



real part of \ where \ = re^^ (dimensionless) 
r for conventional method 
time (seconds) 

real part of \ where \ - X + IZ (dimensionless) 

X for conventional method 
displacement (in) 

imaginary part of \ where \ - X I Z (dimensionless) 

2 for conventional method 

characteristic complex root in the solution of the exact 
equation (d imens ionles s ) 

imaginary part of CL 

real part of CL 

finite difference interval (space or time increment) 

finite difference step in time (sec/time interval) 

damping ratio 
2 

arc sin — dimensionless frequency factor 



^ for conventional method 
^ for MOD 1 method 
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for MOD 2 method 



e. $ 



K 


—In 


r dimensionless 


attenuation factor 


^c 


K 


for conventional 


method 


^1 


K 


for MOD 1 method 




K2 


K 


for MOD 2 method 





K 

\ 



characteristic complex root in the solution to the finite 
difference equations (dimensionless) 

X for conventional method 

X for MOD 1 method 

X for MOD 2 method 

B. 



relative error in 0 



a. At 

jj, for conventional method 

fJL for MOD 1 method 

fJi for MOD 2 method 

K 



relative error in K 



U) 

Co 



1/ for conventional method 
1/ fox MOD 1 method 
J/for MOD 2 method 

natural frequency of vibration (rad/sec) 

CilAf , sampling resolution (rad/time interval) 
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CHAPTER I 



INTRODUCTION 

In the quest for accurate numerical solutions to differential 
equations, considerable effort has been devoted to the determination 
of the errors associated with various finite difference schemes. This 
work has gained importance with the recent widespread use of the 
digital computer. Due to the high value placed on computer time and 
to the limitations on core storage space within the computer, the 
accuracy and efficiency of the numerical method employed are of the 
utmost importance . 

In general, many engineering problems can be classified as 

2 

either boundary value problems or initial value problems. For the 
numerical solution to boundary value problems, the most common method 
used is a form of the central finite difference operators. Several 
forms of these operators can be found in Table III of the Appendix to 
reference 1. Some of these are frequently used to solve engineering 
problems while others seem either to be impractical or to appear 
seldom in current literature. The Hermit ian equations, which are 
included in this table, fall into the last category. For initial 
value problems, numerical methods of solution such as Adams, Runge- 
Kutta, and other forms of the predictor-corrector methods are usually 
employed. The Hermit ian method is also applicable to initial value 
problems but is seldom used. 

The purpose of this paper is to evaluate the merits of the 
Hermitian method when applied to a second order, ordinary, linear 
differential equation. This evaluation will be made by comparing the 
solution to the exact differential equation with the solution to two 
Hermitian finite difference approximations and with the solution to 
the first order central finite difference equation. In this manner, 
the difference in accuracy between these finite difference methods can 
be determined, and the merits of the Hermitian method established. 
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The equation chosen for this analysis is that of the single 

degree of freedom, damped oscillator « This equation was chosen 

because it contains both the first and second derivatives and because 

4 5 

its exact solution is well known*, This appears to be the first 
application of the Hermitian equations to differential equations con- 
taining more than one degree of derivatives,, This problem is generally 
considered an initial value problemo However, the results of this 
analysis are also applicable to boundary value problems since a boundary 
value can be solved using methods of solution normally employed in 
initial value problems,. Also, the oscillator may become a boundary 
value problem if the displacement or velocity is specified at a time 
other than zero. Thus, the choice of this example allows the conclu- 
sions to be applied to both initial value and boundary value problems. 

The author wishes to acknowledge a special thanks to 
Professor Robert E. Ball of the Aeronautics Department of the Naval 
Postgraduate School. The guidance and consultation of Dr. Ball, as 
thesis advisor, were instrumental in the development of this paper. 
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CHAPTER II 



DISCUSSION OF THE CENTRAL FINITE 
DIFFERENCE METHOD AND ITS ASSOCIATED ERRORS 

There are three types of error associated with the finite 

3 

difference solution to engineering problems. One type is called 
truncation or discretization error. The conversion of a differential 
equation to an algebraic equation is, in essence, the approximation 
of a continuous system by a discrete mode. The error involved in 
this conversion is called the truncation error. A graphic example of 
how this error is introduced can be seen by examining the derivation 
of the first order central finite difference equations for the first 
and second derivatives presented in Appendix A. 

In the derivation of the finite difference expression for the 
first derivative, only the first three terms of the Taylor series are 
used. The first four terms are used in the expression for the second 
derivative. The truncation error can be approximated by the first 
unused term in the Taylor series. Therefore, a decrease in truncation 
error can be accomplished by using more terms of the Taylor series in 
the finite difference equation. 

Another type of error is called round-off error. This error 
is generated within the computer and is a function of the amount of 
calculation required and the condition of the numbers involved. Each 
arithmetic computation, such as addition, multiplication, and expo- 
nentiation, has a round-off error associated with it. The computer 
is limited to carrying a specified number of digits. Any amount over 
this number must be dropped for each calculation. To decrease round- 
off error, it is essential to minimize the number of arithmetic calcu- 
lations within the program and to avoid subtraction of large numbers, etc. 

The third type of error is called inherent error. This error, 

unlike truncation error and round-off error, is not caused by the use 

3 

of numerical methods. Inherent error is introduced into the program 
within the data. If the data are not exact, which is usually the case, 
the solution cannot be expected to be exact. 
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The object of any numerical analysis is to obtain as accurate 
an answer as possible without significantly increasing the computation 
time. There are a number of higher order finite difference methods 
available that have less truncation error. Hence^ they provide more 
accurate solutions than the first order central finite difference 
method. Among these is the Hermitian method. These methods vary by 
the number of stations used, the location of the stations, etc. 

When central finite difference methods are applied to an 
engineering problem, the problem is reduced to that of solving a 
matrix. This matrix will usually be a band matrix centered about the 
principal diagonal. The width of the band is dependent upon the number 
of stations included in the finite difference equation. Using the first 
order difference equations listed in Appendix A, the band width would 
be three, and the matrix is called a tridiagonal matrix. If higher 
order difference equations that require more stations are used, the 
band width widens. By widening the band of the matrix, more computa* 
tions may be required than for the tridiagonal matrix. This increases 
the round-off error and the computation time. Furthermore, the tri- 
diagonal matrix will require less core storage than the wider band 
matrix for the same element size. 

Another approach used to decrease the magnitude of the trun- 
cation error involves decreasing the interval between the stations. 

This is equivalent to increasing the degree of the matrix and, there- 
fore, introduces greater round-off error. Other effects of this pro- 
cedure are the increase in computer storage space required and the 
increase in computer running time. 

The Hermitian method considered in this thesis results in a 
tridiagonal matrix but also decreases the truncation error. Conse- 
quently, the gain in accuracy due to the decrease in truncation error 
is not offset by an increase in round-off error. Furthermore, this 
method can be used to advantage on those problems where core storage 
requirements are important. 
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CHAPTER III 



APPLICATION OF THE HERMITIAN METHOD 
TO THE DIFFERENTIAL EQUATION 

The first order central finite difference equations for the 
first and second derivative, as derived in Appendix A, are 

t-i ) 

The Hermitian equations for the first and second derivatives, as 
derived in Appendix B, are 

( ^ I (3.2.) 

where (') denotes the derivative with respect to time. 

As shown in Appendix B, the Hermitian equations have less trun- 
cation error than the first order equations. Note, that in order to 
get this reduction in trucation error it was not necessary to use 
more stations. 

Consider the single degree of freedom, damped oscillator shown 

4 5 

in Figure 1. The differential equation of the oscillator is * 

my"t cyU ky — q (3.3) 

where m is the mass of the system, c is the viscous damping coeffi- 
cient, k is the spring constant, and y is the displacement. Sub- 
stituting the conventional first order central difference equations, 
Eqs. 3.1a and 3.1b, for the first and second derivatives into Eq. 3.3 
leads to 

where ^ has been replaced with^ f to denote the time interval. 
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FIGURE I 




THE SINGLE DEGREE OF FREEDOM, 
DAMPED OSCILLATOR 
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Since the Hermitian equations are not of the same form as 
the conventional equations, a different method of approximation to 
Eq. 3.3 is necessary. Examination of the right side of Eqs. 3.2a 
and 3.2b reveals that Eq. 3.2a contains the first derivative at the 
three stations j+1, j, and j-1. Likewise, the right side of Eq. 3.2b 
contains the second derivative at the same three stations. By adding 
a third equation expressing y at j+1, j, and j~l and modifying one of 
the Hermitian equations, it is possible to solve the differential 
equations at these three stations simultaneously. Equations 3.5, 

3.6, and 3.7 show this procedure. 
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(3.7) 



Adding Eqs. 3.5, 3.6, and 3.7 gives 
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10 m // 4 c ' 10 k 

12 Kj 12 12 V 



(3.8) 
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Note that the right side of Eq. 3*8 will equal zero if ^ 

added to both sides • The conventional approximation for y\ gives 



4^f 






— w ^ 
/ 2 V 



(3.9) 



If Eq. 3.9 is added to Eq. 3.8, the result is 



12 



yj-^i ^fOyj^yj-i 



O 



(3.10) 

This is one form of the Hermitian approximation to Eq. 3.3. 

To get the second form of the approximation, it is necessary 
to write Eq. 3.7 as 






When Eqs. 3.5, 3.6, and 3.11 are added, then 



2m 



Hj + 



2At 



(3.11) 



-yj-i 









(3.12) 



Note that the right side of Eq. 3.12 will equal zero if ™y''j i® 
subtracted from both sides. From the conventional approximation for 

..I I 



m 





(3.13) 
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If Eq, 3.13 is subtracted from Eq. 3.12, the result is 



C 

2Af 






Hj^i Hj-i 



6 \ V V'/ (3.14) 

This is the second form of the Hermitian approximation to Eq. 3.3. 

In this thesis, Eq. 3.10 will be called MOD 2, and Eq. 3.14 will be 
called MOD 1. 

By following the derivation of these two Hermitian approxi- 
mations, it should be clear that a first order central finite 
difference equation for either y’^ or y'V, was used in each deriva- 
tion. For example, the first order central finite difference equa- 
tion for yV, Eq. 3.9, was introduced in the derivation of Eq. 3.10, 
the MOD 2 approximation. This is equivalent to modifying Eq. 3.6 
to be 

2Af 1 (3.15) 

Appendix C shows that the truncation error in Eq. 3.15 is the same 
as the truncation error in the first order central finite difference 
equation for the first derivative, Eq. 3.9. Thus, when the Hermitian 
equations are applied to the differential equation containing more 
than one degree of derivative, accuracy can be gained for only one 
specified derivative at a time. The MOD 1 approximation has the 
accuracy gain in the first derivative, and the MOD 2 approximation 
has the accuracy gain in the second derivative. 

Examination of Eqs. 3.4, 3.10, and 3.14 reveals that the 
only difference in these three equations is in the third term. 
Equations 3.10 and 3.14 allow values for Fj while Eq. 3.4 

does not. This ability to let y vary over the interval j+1 to j-1 is 
the reason truncation error is decreased. If the variation of y 
happens to be linear over the three successive stations j+1, j, and 
j-1, Eqs. 3.10 and 3.14 then both reduce to the conventional equation, 
Eq. 3.4. 
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CHAPTER IV 



COMPARISON OF THE SOLUTIONS 



The three finite difference approximations to the differen- 
tial equation, Eq. 3.3, are repeated here for convenience. 
Conventional 



m 



+ ^ 0 



+ -Hj-f 



(3.4) 



MOD 2 



m 






yj-H -2yj )+ 






MOD 1 



_k_ 

12 



yjH^'%^yj-> ) = ° 
^yj-> )+ 

yj.,^^yj^yj-,\^ ^ 



(3.10) 



^ 0-14) 

The exact solution to the differential equation is of the 

form 

y = A e ^ ^ (4.1) 

where A is an arbitrary constant. Substituting Eq. 4.1 into the 
differential equation, Eq. 3.3, leads to 

/T)a^+ ca ^ ^ =0 

If Eq. 4.1 is to be valid for all values of t, then 



Thus, 



OL = 



m a^-h ca + k = 0 



-c 1 '\/c2 ^4mk 
2 m 



(4.2) 



20 



If 



U) — ■ 



m 



and 






2moo 



, then Eq. 4.2 becomes 



a=— co^ ± /o) V / - ^ 



(4.3) 



^ is the nondimens ional damping coefficient. A value of 
^ less than 1.0 means the system is underdamped, and the motion is 
oscillatory. A value of ^ greater than 1.0 means the system is 
overdamped, and the motion is non-oscillatory. When ^ is exactly 
equal to 1.0, the system is critically damped. For this analysis, 
the special case of critical damping is not considered. 

2 

The solution assumed for the finite difference equations is 

J 



yj = A \ 



(4.4) 



where j is the station number. Solutions for that satisfy 
Eqs. 3.4, 3,10, and 3.14 are derived in Appendix D. For the first 
order finite difference equation, 






4 



60 



(4.5) 

where CO = COA f - 

A discussion of ^ is in order here. The quantity Ci) 
the natural frequency of the oscillator in radians per second. 

Therefore, CO At ^ units of radians per time interval and 

is a measure of the resolution of the solution. A small value of 
CO would mean that the system is being sampled many times per cycle. 
For example, a value ot CO equal to 0.5 would mean the system was 
being sampled every 0.5 radians, or over 12 times per cycle. A value of 
CO over 1.0 would mean that the solution has a very poor resolution, 
and the truncation error would be high. 

Equation 4.1 can be written in the form 



or 



yj = A ( JAt) 






a Af 



J 



(4.6) 
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Comparing Eqs* 4.4 and 4.6 reveals that 



aAt 



vs 



X 



is a proper comparison. However, Eqs. 4.3 and 4.5 show that both 
OL and \ are complex numbers for the underdamped system. There- 
fore. 



C 

where 



aA^ - a. At ± } a At 
— e ' e 



(4.7) 



a^--cjC , O.J = 

and for the first order difference scheme 

X = X ± i z 



where 






The solution \ can also be given in the form 



X nr 
C 



\ = r e' ^ 



where 






and ,9 =/^PC SIN — 



A further simplification can be made by replacing r with 



(? ^ . Thus , 



K = - In r 

and becomes + i f) 

X — ^ — I [y 

r C? e 

Examining Eq. 4.7 and 4.8 reveals that 



(4.8) 



aAt 



vs 



and 



- K 

e 



a. At vs 

I 

are proper comparisons. Let the relative error in both K and Q 
be defined by 



and 



1 / = 



■K 



a- At 



o^At 



Applying the value of Q. and CL ^ to |Lt and 1/ , leads to 

9 ' . K 



fj-- - 



CO 






and 



1 / = 






(4.9) 
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The relative error terms, jJi and J/ , provide an indication 
of the accuracy in the numerical solution. This can be seen by com- 
paring the complete solutions. Substituting Eq. 4.3 into Eq. 4.1 
yields the exact solution 



y = A 






(4.10) 

since f — , f . Substituting Eq. 4.8 into Eq. 4.4 yields the solu- 

tions for the finite difference approximations 



y = A 



[-K ± i e] J 



(4.11) 

Applying the relative error terms jLf and )/ given by Eq. 4.9 to 
Eq. 4.11, the finite difference solutions take the form 



y = A c 



± / /LICOV/ - 



(4.12) 

Thus, a comparison of Eq. 4.12 with 4.10 shows that if jj_ and~l/ 
are both equal to 1.0, Eq. 4.12 is the exact solution. 

Examining Eq. 4.12 reveals that jJi and J/ have a specific 
role in the error analysis. In Eq. A.12 , a measure of the 
rate of attenuation in the finite difference solution. Likewise, 
a measure of the frequency of the oscillations in the 
finite difference solution. Therefore, jy will be an indication of 
the error in attenuation, and jj, will be an indication of the error 
in frequency of the numerical solution. 



I 



!9 



— K 

For the overdamped problem (? does not change. However, 



becomes an attenuating factor and Eq. _4.12 becomes 

y - A 



an attenuating ractor ana hqo D€ 



where 



is given in Eq. 4,9, and U, becomes 
/ 0 




(4,13) 
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The expressions for \ ^ ^ ^ ^ and /Q for Eqs, 3o4g 3,10, and 

3.14 are derived in Appendix D. Table I lists the values of \ , 

K , and ^ , and Table II lists the values for / Q for all three 
numerical approximations. The subscripts 1, and 2 refer to conven- 
tional, MOD 1, and MOD 2 respectively. 

In conclusion, a proper comparison of the solutions to the 
three finite difference approximations with the exaet solution has 
been shown to take the form 



U = 



K 



for the underdamped problem and 



and 



M = r; 



OJ 



K 



e_ 

_e 

thatrmakes the 



o _ and ^ 

L, u OJ 

for the overdamped problem. A computer program 
necessary comparisons was written in the FORTRAN IV language for the 
IBM 360 digital computer. The results obtained from this program 
follow. 
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TABLE I 



LIST OF THE UNDERDAMPED SOLUTIONS 



Xc = 



= 




A PC SIN 1 



-2 




CO - 


^co ~ 7-C0 \ 


/ 





JL 

2. 



X, = 



(/ - -¥]i ;q-\J I - 



( /- ^ ) 



(/ 1 

~[ 7 TT^T^Y~ 



Kj^ zi - ln\ 



6 , = 



X , 



K. zz t In 



ARC SIN ^ 



—2 

cu - ^ 






[\ --^ i G> -\ll - C- 



± 



(/- If] 



( / ^ 6 ^ /f4 ^ "“ Ci^l 






CO 



/A 



0^—arc sin ± 



~2 r^—z / — 

- ; OJ - OJ 



CO 



/ ^ ^ (jS" 
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TABLE II 



LIST OF THE OVERDAMPED SOLUTIONS 



jQ^ = ARC SINH 1 



i — ARC SINH ^ 
i 9^- ARC SINH 1 




I 



Cl) 







~X) 


-2 


/ \ 


U) C 


y “ ^ 


+ ,2 CO 


7 - - 






/ — 


-z 


, 


OJ ( 


0 - 


+ 5 CO 
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CHAPTER V 



RESULTS AMD CONCLUSIONS 



Results 

Figure 2 is a plot of the data in Table III. For this plot, 
^ is equal to 0.0, and the system is undamped. The undamped solu- 
tion will have no attenuation and therefore, the term is meaning- 
less. 

With A = 0,0, the differential equation becomes 



m y'-h k y q 



(5.1) 

Recall that the MOD 2 method was derived by maintaining the Hermitian 
form of the second derivative. As a result, the MOD 2 method pro- 
duces very little error for all values of 60 . Since the first de- 
rivative is not present in Eq, 5.1, the MOD 1 method can be expected 
to have more error than MOD 2, However, it is found that the error 
in MOD 1 is approximately the same as the error associated with the 
conventional method. 

The significant gain in accuracy by the employment of MOD 2 
means a larger time interval could be used to approximate the system 
without loss of accuracy. For example, if the natural frequency of 
the system equaled 180 radians per second, a time interval of 0.005 
seconds would describe the system with 0,2 per cent error in 
using MOD 2. To achieve this accuracy with the conventional method, 
a time interval of 0.0015 seconds must be used. This means that the 
time interval can be increased by a factor of three by using the 
MOD 2 approximation for this problem. 

As the damping is increased, the coefficient of the first 
derivative increases, and the first derivative begins to have a sig- 
nificant influence upon the response. The MOD 1 method should become 
more accurate as ^ increases since it is derived by maintaining the 
Hermitian form of the first derivative. Figures 3 and 4 are plots of 
the data presented in Table IV with ^ equal to 0,9, Since ^ is 
less than 1,0, the underdamped comparisons are used. As anticipated. 
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a»At RAD/ TIME INTERVAL 

FIGURE 2 ^ vs SAMPLING RESOLUTION FOR DAMPING = 0.0 



TABLE III 

RELATIVE ERROR FOR K ® 0,0 



CO 


H-c 






0,1 


1.00042 


1 

0,99958 


1.00000 


0.2 


1.00167 


0.99834 


1.00000 


0.3 


1.00379 


0.99629 


1.00002 


0,4 


1.00679 


0.99345 


1,00005 


0,5 


1.01072 


0,98987 


1.00013 


0,6 


1,01564 


0.98559 


1.00027 


0.7 


1.02163 


0.98065 


1.00051 


0.8 


1,02879 


0.97513 


1.00088 


0.9 


1,03726 


0.96907 


1.00142 


1.0 


1,04720 0.96255 ^ 


- 

1 00219 



29 



MOD 



--CO 




- • lO 



30 



cuAt RAD /TIME INTERVAL 

FIGURE 3 /I vs SAMPLING RESOLUTION FOR DAMPING = 0.9 



fsi 
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FIGURE 4 V vs SAMPLING RESOLUTION FOR DAMPING = 0.9 



TABLE I¥ 



RELATIVE ERROR FOR ^ * 0.9 



CJ 








h 






0.1 


1.0033 


0.99776 


1.0022 


1.00051 


1.0022 


0.99914 


0.2 


1.0116 


0.99072 


1.0044 


1.00203 


1.0078 


0.99643 


0.3 


1.0248 


' 

0.97780 


1.0096 


1.00456 


1.0174 


0.99151 


0.4 


1.0466 


0.95681 ! 


1 1.0169 


1.00809 


1.0319 


0.98364 


0.5 


1.0771 


0.92353 1 


1 1.0273 


1.01260 


1.0515 


0.97153 


0.6 


1.1190 


0.86943 I 


1 1.0409 


1.01804 


1.0781 


0.95291 


0.7 


1.1770 


0.77470 \ 


i 1.0571 


1.02432 ! 


il.1133 


0.92357 


0.8 


1.2605 


0.57647 i 1.0785 


1.03123 i 


1 1.1608 


0.87490 


0.9 


1.3913 


1 

0.0 1.1044 


1.03841 


1.2259 


0.78621 


1.0 


1.6360 


1 

0.0 1 1.1375 


1.04507 


1.3230 


0.58758 


} 
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the MOD 1 approximation is seen to be the most accurate. For example, 
in Figure 3 at a resolution of 0.8 the MOD 1 approximation has 3 per 
cent error in ^ • For this same resolution, the conventional approxi- 
mation has 42 per cent relative in Q . 

Figures 5 and 6 are plots of the relative error versus ^ 
for a constant resolution of 0.5 as listed in Table V. Figure 5 
shows that one of the two Hermitian approximations always has a better 
value of jJL than the conventional approximation. For low values 
of ^ the MOD 2 approximation is best, and for high values of ^ the 
MOD 1 approximation is best. This means that throughout the range 
of ^ , greater accuracy in the frequency can be obtained by using one 

of the Hermitian approximations. 

The error in attenuation, J/ , is shown in Figure 6. Note that 
for the lightly damped system, where the attenuation is small, the 
conventional approximation has the least error. However, this does not 
imply that the total solution given by Eq. 4.12 is most accurate when 
the conventional approximation is employed since Figure 5 shows that 
the conventional approximation has the largest error in jLi for the 
lightly damped system. For values of ^ between 0.45 and 0.65, the 
MOD 2 approximation is best, and for values of ^ greater than 0.65, 
the MOD 1 approximation is best. Thus, Figure 6 shows that the Her- 
mitian approximations have less error in attenuation for values of ^ 
greater than 0.45. 

Analysis was not made for the critically damped oscillator. 
However, as ^ approaches 1.0 from either side, Figure 5 shows that 
jJi becomes very large. In order to understand what happens to the 
numerical solution for ^<Cl.l, it is necessary to make other 

comparisons. For example, the exact solution, Eq. 4.10, and the 
solution for the finite difference approximations, Eq. 4.12, can both 
be given in the following form 

r-^ £ ; V / - {;■“] CO Af ^ 



and 



y- A c 
y-A e 



^ j/ t /'jU-V/ - ^ ^ 
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FIGURE 6 V vs DAMPING FOR SAMPLING RESOLUTION 



TABLE ¥ 

RELATIVE ERROR FOR W = 0.5 
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1.01072 




0.98987 




1.00013 ! 
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1.000 


1.01125 


.9600 
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.9780 


1.00074 1 


.2 


1.004 


1.01282 


.9630 


0.99259 


.9830 


1.00255 1 


.3 


1.008 


1.01531 


.9666 


0.99598 


.9873 


1.00551 1 
1.00948 


.4 


1.0140 . 


1.01849 


.9725 


1.00070 


.9920 


.5 


1.022 


1.02184 


.9792 


1.00665 


1.0004 


1.01418 


.6 


1.0317 


1.02419 ! 


.9880 


1.01356 ■ 


1.0096 


1.01890 


.7 ! 


1.0440 


1.02254 
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1.02067 


1.0211 


1.02181 
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1.0630 


1.00712 


1.0120 


1.02527 


1.0350 


1.01692 


.9 


1.0771 
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1.0273 


1.01260 , 


1.0515 


'0.97153 


1,1 


1.1244 


1.34268 j 


1.0680^ 


1,17106 


1.0950 


1.25556 


1.2 


1.1552 


1.29560 1 


i 1.0940 


1.17212 


1.1257 


1.23259 
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1.1926 
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1.19704 j 
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1.34629 
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1.40661 


1.2101 
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1.2517 J 
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for ^ Thus a plot of,^ I ^ and|^\y/«- ^^ versus ^ shows 
the error in the oscillatory part of the solution for the underdamped 
system. For the overdamped system, ^ ^ loO, Eqs. 4.10 and 4.12 



become 



= 4 J 



and 



= 4 c[-^^ ' J 






Thus, for this range ^ ^ plot of ^ | and ^ J versus 

will be the equivalent comparison for the overdamped system. 

Examination of the equations in Table I reveals that the 
numerical solution will be critically damped for the values of ^slightly 
less than 1.0. Therefore, to show the critical damping point for each 
finite difference approximation it is necessary to use the comparisons 
for the overdamped system for values of .95^^Cl.O. Equation 4.13 
defines jJi for the pverdamped system as 



therefore. 









/& 



I 



CO 



(5.2) 

Using Eq. 5.2, it is possible to evaluate the overdamped comparisons 
for values of ^ <C 1.0. 

Figure 7 is a plot ofjJi\Jl - ^^or jJiyJ j for values of ^ 
between 0.9 and 1.1 with CO * 0.5. The MOD 1 approximation has its 
critical damping point closest to 1.0, and the conventional approxima- 
tion has its critical damping point furthest from 1.0. As CO is de- 
creased and the system is sampled more often, the critical damping 
points of all the numerical solutions will approach 1.0. However, the 
Hermitian approximations will always be critically damped closer to 
1.0 than the conventional approximation. The fact that' each numerical 
solution has its own critical damping point that is less than 1.0 
explains the large values of jJi in this region. 

For the overdamped system a problem arises when ^ (X) is greater 
than 1.0. Examination of the equations in Tables I and II reveals an 
unstable zone for these finite difference equations at this point. 
Although ^ Cl) = 1«0 is essentially the limit for numerical stability. 
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FIGURE 7 INVESTIGATION OF THE REGION NEAR CRITICAL DAMPING 



the accuracy of the difference approximations becomes very poor 
before the unstable zone is reached. As the stability limit is 
approached 9 both jJi and J/ approach infinity. Figure 8 shows the 
stability limits on all three finite difference approximations. This 
plot was made using the equations in Tables I and II . 

Conclusions 

From the results presented hercj, the conclusion can be made 
that the Hermitian equations can be used to obtain more accurate 
numerical solutions for the same size increment than can be obtained 
by first order central finite difference equations when applied to 
differential equations containing derivatives of more than one degree. 
Further, this Increase in accuracy is not accompanied either by an 
increase in complexity of the resulting matrix when applied to boundary 
value problems or by a significant increase irjL computation time when 
applied to initial value problems. Consequently, use of the Hermitian 
equations would be most advantageous when computer core storage space 
is limited and accurate results are required. 
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FIGURE 8 DAMPING vs SAMPLING RESOLUTION NEAR STABILITY LIMIT 
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APPENDIX A 



THE DERIVATION OF THE CONVENTIONAL EQUATIONS 
FOR THE FIRST AND SECOND DERIVATIVE 



The first order central finite difference equations can be de- 
rived from the Taylor series. The stations are lableled station 
j+1, station j, and station j-1. The interval between stations j+1 
and j is called 

The Taylor series for V is 

1 , ^ 

II A / A A A ^ 



2- 

J± U f 

S! V 



For Yj the Taylor series becomes 

yj- = ^ - ^ rf y/- 



^ ” 



—y, + 

S/ V 



A^ ✓// 



Subtracting Eq. A.l from Eq. A. 2 gives 



(A.l) 



(A. 2) 



/ A3 A^" ^ 

yj.,- yj-i = ^^yj ^ t^j ^ "" 



iJ 



U 



(A. 3) 



Solving for y'^ leads to 






_ / 



y,.. - y 



A 2 /// A f JET 

^y,- y,^- - 



(A. 4) 



This is the first order central finite difference equation for the 
first derivative. 

Adding Eqs. A.l and A. 2 gives 




/ 

A^ 




-l^y 

A^ 

36o^j 




A1 



(A, 5) 



Equation «.5 is the first order central finite difference equation 
for the second derivative. 
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The truncation error in each equation is usually approximated by 
the first term in the Taylor series that is not used in the difference 
approximation. For example, the error in Eq. A. 4 is 

A ^ 

- 6 - 



and the error in Eq, A. 5 is 



m 



tl 
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APPENDIX B 



THE DERIVATION OF THE HERMITIAN EQUATIONS 
FOR THE FIRST AND SECOND DERIVATIVES 



Taking the derivative of Eq. A.l in Appendix A gives 
/ 3 j]? 






/// 

■y, ^ 



A 



1 ' 6 V 



-h 



A 






^4 II 

Likewise Eq<, Ao2 becomes 



y, + 



Hj-> 



yj- ^ yj^ 

^ 7.4 yj 



Adding Eqs. B.l and B.2 gives 



A^ 0 

— U + * » • 

/XO^J 



3 



A 2 /// ^ 

^yj~ -tHj 

/xoyj ^ 



(B.l) 



(B.2) 



^ -z 



yjH ^ yj-' ^^y' ^ ^y< ' 



t/ 



A /// A 



Dividing Eq. B.3 by 6 and solving for A ^ y' ' ' leads to 

6 

A^u 44u' - a y ' ^ y' 1- A’L - 

6 V V-J-) V 

From Eq. A. 4 of Appendix A, 



= -^(c 






/// 



'yj^ztvj^r y/-.i-^/- 



A 



L /- 1 1 6 ^ ^ 



Thus, substituting Eq. B.4 into Eq. A. 4 gives 



yJ = 4 



yj..- yj-}-i[yjH-^yy^ V-l 



A^ . , 



(B.3) 



(B.4) 



(A. 4) 



(B.5) 
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Solving for yj^.j " i gives 



= jfe "S' ' 



s-.; 









lo 


1 

t 


Equation 


B.6 is 


the 


Hermitian 


equation 


The 


second 


derivative of 


equation . 


// 


// 


* . . "! 


A^,i^ 








A <dj f 










^ 0 


A ^ HI 








z^yj ^ 


1X0 yj 


Likewise 


Eq • A. 2 


becomes 




// 


ft 


/// 


132 








A 9j- 


X^- 








A’yZl 

XfV 


IXO^J 


Adding Eqs. B.7 


and 


B.8 gives 




// 


// 




// 




u 

Vi-t 


" yj-> 


= + 


^^yj 



(B.6) 












" /1-yj 

Dividing Eq. B.9 by 12 and solving for ^ ^ y yields 

12 



A-2 / / // If 



From Eq. A. 5 of Appendix A, 



// 












u" = J- 



^jH -^yj^ 

J6<? 

Substituting Eq. B.IO into Eq. A. 5 gives 



O-^Tryj 



(B.7) 



(B.8) 



(B.9) 



(B.IO) 



(A. 5) 



-L 

V A^ 



/ 



// 



" ^^-'1 

■h 



/\^ ^ 

^ y -h * * * 






(B.ll) 
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Then, solving for (y 



j+i ■ <■ 5’j-P 



y. 


-^y. 


t y 


II 


I '' 

\y, 


// 

■h to y -h 


y, 




j 




1 U 1 




J 


j-‘\ 



A ^ SL 

u f- * ^ * 



(Ba2) 



Equation B.12 is the Hermitian equation for the second derivative. 

A comparison of the conventional equations in Appendix A with 
the Hermitian equations derived here reveals that the Hermitian 
equations have less truncation error. This is to be expected since 
the error term in each first order equation was used to derive the 
respective Hermitian equation. 
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APPENDIX C 



DISCUSSION OF THE TRUNCATION ERROR 
IN THE MODIFIED HERMITIAN EQUATION 



In the MOD 2 method of analysis, the Hermitlan equation for 
the first derivative is modified. Equation B.6 is the Hermitlan 
equation for the first derivative 









This equation can be written 



A ^ 

— -= — u ^ 

3iO V 

The conventional equation for y'^ can be written 



(B.6) 



(C.l) 



C 






i/-' 



c ' A-2 

= T'^j *-n:^j 



Addition of equations C.l and C.2 gives 



(C.2) 



C 



^A 



- 






_ c 






J+! 

Al 

/X 



/// A ^ 

^ ~ 360^J 



(C.3) 

Equation C.3 is the modified Hermitlan equation for the first 
derivative that was used in the MOD 2 approximation. The truncation 
error in this equation is 



/X 






/// 



This truncation error is one-half the truncation error in the first 
order equation for the first derivative. Therefore} in the MOD 2 
approximation to the differential equation, there is very little gain 
in accuracy associated with the first derivative. 
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Likewise, for the MOD 1 approximation, the Hermitian equation 
for the second derivative must be modified The modified equation is 



/m 

Ai 





i- 










IX V IfO 



The truncation error, 



(C.4) 



A^ 

1 

is equal to the truncation error of the first order difference 
equation for the second derivative. Therefore, in the MOD 1 approxi- 
mation, there is no gain in accuracy associated with the second 
derivative. 




48 



APPENDIX D 



SOLUTIONS FOR THE 
FINITE DIFFERENCE APPROXIMATIONS 



The conventional approximation will be used as an example for 
this appendix. The conventional finite difference expression for 
the differential equation is 



* Kj] ^ 2Kt\ 



*kyj = 0 



(D.l) 



When the assiomed solution 



yj = A \ 



j 



is substituted into Eq. D.l, the result is 

r ^ 



UAtT 






= 0 

For an arbitrary solution. 



m 

(A 



Xc ^x 



-h 

^ I 2C^f 



+ /< X^ =0 



or 



r m c 
UAtf 






2 m 

2fl 



K 



m 



ZAt_ 



- 0 




By using the quadratic equation. 



1 1 

■•^ 1 ^ ■ 


Ml 


-£i;h 








-f- 



C 

2 Af 



] 
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or 






[z/77 j 



2 



m 



Af 

2 



When the substitutions 



GJ - 




; C~ 

are made, becomes 



and CO-U)Af 



X.- 



i oa/ C - 






[/ ■*■ Caj j 

To make ^ ^ complex, take -^y outside the radical. Then 



Xc- 



^ ± / ojA/ ' - 



[/ ^ ^ c ] 

The solution X Is now of the form 

c 

X^= X + / Z 
To get X Into the form 

■ x.= . . 

the following relationships are used. 



(D.2) 



=: V?7 

Therefore, 



and 






0^ = ARC SIN 

_2f 



or 




Similarly, 



^c- “ I A ~ ^ ^ 

^ I + 



(D.3) 



0^— ARC SIN 1 ^ - -Z 

^ 






(D.4) 



For the overdamped system, i ^ must be found. From Eq. D.4, 




/ SIN Q - SINH I Q 

Therefore, 




-ARC 



SI NH t 



- d)-^ 

V / - 



(D.5) 

The equations for \ , K , ^ , and i ^ for MOD 1 and MOD 2 
are derived in a similar manner. 
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